1. Model eigenvalues and impulse response

Eigenvalues play a crucial role in the dynamics of models, encapsulating the adjustment process within a model. Their sign and magnitude dictate whether a model will stabilize, diverge after a disturbance, or if such a disturbance will induce oscillations.

For macromodels, which are systems of (non)linear difference equations, the initial step involves linearizing the model by computing Jacobi matrices. Since the model may have lags longer than 1 so the linear model can be of an order higher than one, it is transformed into a first-order system by creating a companion matrix. The stability, marginal stability, or instability of the system is assessed based on the eigenvalues of this companion matrix.

The Modelflow library provides a comprehensive set of tools for calculating a model’s eigenvalues, facilitating the assurance of the model’s dynamic stability and identifying potential sources of instability.

1.1. Load the model

mpak,bline = model.modelload('..\models\pak.pcim', 
             alfa=0.8,run=True,keep= 'Baseline',start=2023, end=2100)
mpak.basedf = bline

file read:  C:\modelflow manual\papers\mfbook\content\models\pak.pcim

1.2. Make some scenarios

To provide insight into how shocks impact the model results each scenario shocks the oilprice in 2025. The impulse response function then illustrates the dynamic response of the model.

for oilshock in [0.01 ,5. ,10. ]: 
    oilshockdf = bline.upd(f'<2025> WLDFCRUDE_PETRO + {oilshock}')
    _ = mpak(oilshockdf,2023,2100,keep=f'${oilshock} increase in oil prices 2025',alfa=0.7) 

1.3. .eigenvalues_plot() Display polar diagram of eigenvalues

This is a method to get an overview of the eigenvalues. The .eigenvalues_plot() method first finds all eigenvalues in current time frame, then plots a polar diagram of the eigenvalues for selected years.

Inspection of the polar diagrams reveals indicates that most of the eigenvalues are less than one, and that those that have imaginary roots (those that appear off of the 0 degree line), only a few are outside of the unit circle, indicating exploding osiliation. However not all years have such values. which could be a source of model uncertainty (these are hard to see in the figure, but are identified below). However the real values outside the unit circle indicates that the model will have exploding features - more on that later.

fig = mpak.eigenvalues_plot(periode=[2023,2024,2099,2100],size = (6,6));

Note

fig can be exportet to file with:

fig.savefile(‘polar.pdf’)

If another format then pdf is needed, .pdf can be changed to ´.svg|.png|.eps|.jpeg´.

Inspection of the polar diagrams reveals indicates for Pakistan model:

• Real eigenvalues larger than 1, suggesting an exploding response to shocks.

• There are instances where the eigenvalues indicate a response that is both exploding and oscillating.

Nevertheless, for the explosion and oscillation effects to be significant, they need to persist over several years, which means the occasional occurrences might not have substantial effects.

1.4. the method .eigenvalues_show() Interactive exploration of eigenvalues and vectors

To enable closer analysis the .eigenvalues_show() function provides a interactive interface to the analysis of models eigenvalues and vectors. It will generate and display a polar plot of eigenvalues and their corresponding eigenvectors for a selected year. The user can interact with the plot via a dropdown for year selection, a slider for eigenvalue selection, and a button to toggle additional plot details. The polar plot dynamically updates to reflect the selected eigenvalue, displaying its magnitude and phase. Additional information about the selected eigenvector is also displayed, facilitating a detailed temporal analysis of the eigenvalues.

by pushing the open plot widget button a widget to look at the impulse response for the variables which contributes with weights in the eigenvector corosponding to the selected eigenvalue.

mpak.eigenvalues_show() 

The plot widget

The impulse response for the one variable in the first Eigenvector

1.5. Impulse response

The impulse response can be visualized with the .keep_plot function. Below, observe the impulse response for selected variables. Initially, we present the response for the maximum projection period to illustrate long-term behavior, followed by a shorter period to highlight oscillations in the short run.

For many variables, the response amplifies over the long run, consistent with real eigenvalues above 1. Similarly, some variables show a diminishing oscillation, corresponding with the presence of eigenvalues with imaginary components.

1.5.1. Long run absolute impulse response functions

mpak.keep_plot(' PAKNYGDPMKTPCN PAKBNCABFUNDCN PAKGGDBTTOTLCN PAKCCEMISCO2TKN',
diff=True,start=2023,end=2100,samefig=1,legend=1,title='Impact long run');

1.5.2. Long run relative impulse response functions

In interpreting this, it’s crucial to acknowledge that the local price level in the baseline grows significantly. Observing the percentage difference from the baseline reveals that measured in this way the response does not explode. So the real impact of the impulse does not explode.

mpak.keep_plot(' PAKNYGDPMKTPCN PAKBNCABFUNDCN PAKGGDBTTOTLCN PAKCCEMISCO2TKN',
diffpct=True,start=2023,end=2100,samefig=1,legend=1,title='Impact full projection period');

1.5.3. Short run

 mpak.keep_plot(' PAKNYGDPMKTPCN PAKBNCABFUNDCN PAKGGDBTTOTLCN PAKCCEMISCO2TKN'
,diff=True,start=2023,end=2040,samefig=1,legend=1,title='Impact short run');

1.6. Summary:

This workbook focuses on analyzing eigenvalues and impulse responses to understand the stability and behavior of a dynamic model. Key aspects include:

• Eigenvalues and Eigenvectors: Eigenvalues and eigenvectors are useful for assessing the system’s stability. Eigenvalues highlight the system’s potential for growth, decay, or oscillation, while eigenvectors determine the direction of these dynamic changes.

• Impulse Response Analysis: This section explores the impulse responses of selected variables, showcasing their behavior through long-term trends and short-term oscillations. This method provides an understanding of how variables respond to external shocks.

• Visualization:The workbook presents tools for the detailed examination of eigenvalues and impulse responses.

• On the Pakistan Model: Initially, the model appears to exhibit an exploding impulse response. However, this primarily reflects price developments. When measured in relative terms, the impulse response does not explode but approaches zero.

1.7. Appendix: Compainion matirix

1.7.1. A model:

A normalized model with:

• \(\textbf n\) number of endogeneous variables

• \(\textbf k\) number of exogeneous variables

• \(\textbf r\) max lag of endogeneous variables

• \(\textbf s\) max lag of exogeneous variables

• \(t\) time frame (year, quarter, day second or another another unit)

can be written:

(1.1)\[\begin{eqnarray} y_t^1 & = & y_t^2...,y_{t}^n...y_{t-r}^1...,y_{t-r}^n,x_t^1...x_{t}^k,...x_{t-s}^1...,x_{t-s}^k) \\ y_t^2 & = & y_t^1...,y_{t}^n...y_{t-r}^1...,y_{t-r}^n,x_t^1...x_{t}^k,...x_{t-s}^1...,x_{t-s}^k) \\ \vdots \\ y_t^n & = & y_t^1...,y_{t}^{n-1}...y_{t-r}^1...,y_{t-r}^n,x_t^1...x_{t}^r,x..._{t-s}^1...,x_{t-s}^k) \end{eqnarray}\]

Written in matrix notation where \(\textbf{y}_t\) and \(\textbf{x}_t\) are vectors of endogenous/exogenous variables for time t

(1.2)\[\begin{eqnarray} \textbf{y}_t & = & \textbf{F}( \textbf{y}_t \cdots \textbf{y}_{t-r},\textbf{x}_t \cdots \textbf{x}_{t-s}) \end{eqnarray}\]

The functions are normalized, meaning:

• Each endogenous variable is on the left hand side one time - and only one time.

• An endogenous variable without lags can not be on the right hand side in an equation, which has the variable on the left hand side.

1.7.2. The derivatives

We can express the matrices of derivatives with respect to the endogenous and exogenous variables like this:

\[\begin{eqnarray*} A & = & \frac{\partial F}{\partial y_t^T} \\ \\ E_i & = & \frac{\partial F}{\partial y_{t-i}^T } \hspace{5 mm} i=1, \cdots , r\\ \\ F_j & = & \frac{\partial F}{\partial x_{t-j} ^T} \hspace{5 mm} j=0, \cdots , s\\ \\ \end{eqnarray*}\]

1.7.3. Constructing the companion matrix

To calculate the effect of small perturbations around a solution the derivative matrices can be used.

For simplicity a system with \(r\) (max lag of endogenous variables) of 2 and \(s\) (max lag of exogenous variables) of 1 we have the linearized model:

\(\Delta \textbf{y}_t = \textbf{A}\Delta \textbf{y}_{t} + \Delta E_1 \textbf{y}_{t-1} + \Delta E_2 \textbf{y}_{t-2} + E_3 \Delta \textbf{y}_{t-3} + F_0 \Delta x_t + F_1 \Delta x_{t-1}\)

Rearranging the equation to isolate \(\Delta \textbf{y}_{t}\) on the left side:

\((\mathbb{I}-\textbf{A}) \Delta \textbf{y}_{t} = E_1 \Delta \textbf{y}_{t-1} + \Delta E_2 \textbf{y}_{t-2} + \Delta E_3 \textbf{y}_{t-3} + F_0 \Delta x_t + F_1 \Delta x_{t-1}\)

\(\Delta \textbf{y}_{t} = (\mathbb{I}-\textbf{A})^{-1} E_1 \Delta \textbf{y}_{t-1} + (\mathbb{I}-\textbf{A})^{-1} \Delta E_2 \textbf{y}_{t-2} + (\mathbb{I}-\textbf{A})^{-1} \Delta E_3 \textbf{y}_{t-3} + (\mathbb{I}-\textbf{A})^{-1} F_0 \Delta x_t + (\mathbb{I}-\textbf{A})^{-1} F_1 \Delta x_{t-1}\)

This is a 3th order system of difference equations. In order to be able to find the gain stability of the system it has to be rewritten as 1st order system of difference equations. That is a system with only one lag. Fortunately there is a standard way to do this. It runs as follow:

Using the state vector: \(\textbf{z}_t = \begin{bmatrix} \textbf{y}_t \\ \textbf{y}_{t-1} \\ \textbf{y}_{t-2} \end{bmatrix}\) and \(\textbf{w}_t = \begin{bmatrix} \textbf{x}_t \\ \textbf{x}_{t-1} \end{bmatrix}\)

The system can be expressed as:

\(\Delta \textbf{z}_{t} = \underbrace{\begin{bmatrix} (\mathbb{I}-\textbf{A})^{-1}E_1 & (\mathbb{I}-\textbf{A})^{-1}E_2 & (\mathbb{I}-\textbf{A})^{-1}E_3 \\ \textbf{I} & \textbf{0} & \textbf{0} \\ \textbf{0} & \textbf{I} & \textbf{0} \end{bmatrix}}_{\textbf{C}} \Delta \textbf{z}_{t-1} + \left[\begin{matrix}(\mathbb{I}-\textbf{A})^{-1} F_{0} & (\mathbb{I}-\textbf{A})^{-1} F_{1} \\{0} & {0} \\{0} & {0} \\\end{matrix}\right] \Delta \textbf{w}_t\)

This system is a first order system - and the stability can be evaluated using the matrix:

\(\textbf{C} = \begin{bmatrix} (\mathbb{I}-\textbf{A})^{-1}E_1 & (\mathbb{I}-\textbf{A})^{-1}E_2 & (\mathbb{I}-\textbf{A})^{-1}E_3 \\ \textbf{I} & \textbf{0} & \textbf{0} \\ \textbf{0} & \textbf{I} & \textbf{0} \end{bmatrix}\)

also known as the companion matrix

If the systematic in construction the companion matrix with max lag endogenous variable of 2 can be extended to any numbers of lags.

1.7.4. Stability and the eignvalues of the companion matrix

The dynamic behavior can be evaluated by looking at the eigenvalues \(\textbf{e}_{t}\) of the companion matrix \(((I-\bar A)^{-1}\bar E )\). Note that \(\textbf{e}_{t}\) can are complex numbers.

• If all \(\lvert \textbf{e}_{t} \lvert < 1\) the system will converge.

• If at least one of the eigenvalues is larger than one, the system will amplify.

• If at least one \(\textbf{e}_{t}\) has an imaginary part the system will oscillate

  • dampened if all \(\lvert \textbf{e}_{t} \lvert < 1\) or

  • amplifying if one \(\lvert \textbf{e}_{t} \lvert > 1\).

The eigenvalues and associated eigenvectors of a model can be calculated using the method .get_df_eigen_dict(), which returns a dictionary where the key is the year and the value is a dataframe the first row of which is the eigenvalues of the model in that year and the following rows are the eigenvectors.

1.7.5. Companion matrix in modelflow

To construct the companion matrix and calculating the eigenvalues and eigenvectors modelflow goes through a number of steps:

• finding expressions for all partial derivatives with respect to endogeneous variables

• creating a model to calculate the partial derivatives.

• using this model calculate all partial derivatives for all lags

• use the calculated partial derivatives to construct a companion matrix for all year

• calculate the eigenvalues and eigenvectors.

In modelflow these steps are performed in a special class newton_diff- as is was constructed to enable solving the model using the newton method. The model class (from which this notebook has the instance mpak wraps the relevant functions. These are:

• .get_eigenvalues returns a dictionary where the key is the year and the value is a dataframe containing s the eigenvalues of the model in that year.

• .get_eigenvectors returns a dictionary where the key is the year and the value is a dataframe the first row of which is the eigenvalues of the model in that year and the following rows are the eigenvectors.

• .eigenvalues_show An interactive widget to display eigenvalues

• .eigenvalues_plot Returns a polar plot of the eigenvalues in a selected year

• .get_eigen_jackknife_df Compute and cache eigenvalues with each endogenous variable excluded one at a time. Returns a tall dataframe with the results

• .jack_largest_reduction Identifies the largest reduction in eigenvalue magnitude for a specified period and optionally focuses on eigenvalues with imaginary parts

• .jack_largest_reduction_plot Creates an interactive Plotly plot to visualize the reduction in eigenvalue magnitude across different exclusions

• .get_newton Returns an instance of newton_diff class which is used for the eigenvalue calculations.

Using some of these functions may require careful inspection of the help files and an understanding of the theory behind.

When eigenvalues has been calculated once for a model the results are stores in a instance of the class newton_diffclass. The instalse is called .stability_newtonand the calculations are not repeated every time the eigenvalues are used.